Description

I begin this exercise by loading a GeoJSON file from the City of Austin’s Open Data Portal, which shows the locations of 20 community recreation centers throughout the city.

ATX_rec <- st_read("https://data.austintexas.gov/resource/8dff-2vkt.geojson")
## Reading layer `8dff-2vkt' from data source `https://data.austintexas.gov/resource/8dff-2vkt.geojson' using driver `GeoJSON'
## Simple feature collection with 20 features and 10 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -97.81079 ymin: 30.18567 xmax: -97.63703 ymax: 30.43875
## geographic CRS: WGS 84

Next I import the OpenStreetMap data for Austin’s street network.

TxStateCent <- "+proj=lcc +lat_1=30.11666666666667 +lat_2=31.88333333333333 +lat_0=29.66666666666667 +lon_0=-100.3333333333333 +x_0=700000 +y_0=3000000 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs"

austin_street_features <- opq(bbox = 'Austin TX USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_sf()

austin_streets <- austin_street_features$osm_lines %>%
  st_transform(crs = TxStateCent)
ggplot(austin_streets) +
  geom_sf() +
  theme_map()

Open Trip Planner setup

In order to create isochrones for various modes, I first installed the Open Trip Planner Java utility.

After installing the utility, I’ll finally launch OTP and connect to the open utility.

otp_setup(otp = path_otp, dir = path_data, memory =1024)
otpcon <- otp_connect()

Plotting isochrones

Once connected to OTP, I can generate isochrones. I’ve chosen to examine walking and biking access to the community centers in my dataset, and will begin by establishing 5-minute and 10-minute isochrones for these modes.

iso_5min_walk <- 
  otp_isochrone(otpcon = otpcon, fromPlace = ATX_rec, 
                mode = "WALK", cutoffSec = 300) %>%
  st_transform(crs = TxStateCent) %>%
  mutate(mode = "walk")
## Warning in otp_isochrone(otpcon = otpcon, fromPlace = ATX_rec, mode = "WALK", :
## Failed to get isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 300},"id":"fid-1272ace8_17501947ade_-7ff8"}]}
iso_5min_drive <- 
  otp_isochrone(otpcon = otpcon, fromPlace = ATX_rec, 
                mode = "CAR", cutoffSec = 300) %>%
  st_transform(crs = TxStateCent) %>%
  mutate(mode = "drive")

iso_5min_bike <- 
  otp_isochrone(otpcon = otpcon, fromPlace = ATX_rec, 
                mode = "BICYCLE", cutoffSec = 300) %>%
  st_transform(crs = TxStateCent) %>%
  mutate(mode = "bike")
## Warning in otp_isochrone(otpcon = otpcon, fromPlace =
## ATX_rec, mode = "BICYCLE", : Failed to get isochrone
## with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 300},"id":"fid-1272ace8_17501947ade_-7fd0"}]}
iso_10min_walk <- 
  otp_isochrone(otpcon = otpcon, fromPlace = ATX_rec, 
                mode = "WALK", cutoffSec = 600) %>%
  st_transform(crs = TxStateCent) %>%
  mutate(mode = "walk")

iso_10min_drive <- 
  otp_isochrone(otpcon = otpcon, fromPlace = ATX_rec, 
                mode = "CAR", cutoffSec = 600) %>%
  st_transform(crs = TxStateCent) %>%
  mutate(mode = "drive")

iso_10min_bike <- 
  otp_isochrone(otpcon = otpcon, fromPlace = ATX_rec, 
                mode = "BICYCLE", cutoffSec = 600) %>%
  st_transform(crs = TxStateCent) %>%
  mutate(mode = "bike")

otp_stop()
## [1] "SUCCESS: The process \"java.exe\" with PID 33624 has been terminated."

I bind these sets of isochrones to create compound models, which will aid in visualization later. While I created isochrone data for driving, I will choose only to explore walking and biking data.

iso_all_modes <- rbind(iso_5min_drive, iso_5min_bike, iso_5min_walk)

iso10_all_modes <- rbind(iso_10min_drive, iso_10min_bike, iso_10min_walk)

iso10_two_modes <- rbind(iso_10min_bike, iso_10min_walk)

iso_two_modes <- rbind(iso_5min_bike, iso_5min_walk)
right_side <- st_bbox(iso10_two_modes)$xmax
left_side  <- st_bbox(iso10_two_modes)$xmin
top_side <- st_bbox(iso10_two_modes)$ymax
bottom_side <- st_bbox(iso10_two_modes)$ymin

ggplot(iso_two_modes) +
  annotation_map_tile(zoomin = 0, type = "cartolight", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = ATX_rec) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area reachable within 5 minutes",
                       labels = c("By bike", "By foot"),
                       option = "plasma") +
  theme_map() + 
  theme(legend.position = "bottom") +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

right_side <- st_bbox(iso10_two_modes)$xmax
left_side  <- st_bbox(iso10_two_modes)$xmin
top_side <- st_bbox(iso10_two_modes)$ymax
bottom_side <- st_bbox(iso10_two_modes)$ymin

ggplot(iso10_two_modes) +
  annotation_map_tile(zoomin = 0, type = "cartolight", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = ATX_rec) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area reachable within 10 minutes",
                       labels = c("By bike", "By foot"),
                       option = "plasma") +
  theme_map() +
  theme(legend.position = "bottom") +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

A dramatic increase in service area clearly occurs when increasing the isochrone time from 5 minutes to ten minutes.

Comparing isochrone areas

By calculating and comparing the areas of these two sets of isochrones, we can figure out just how much the service areas for the 20 recreation centers increases from a 5-minute-distance to a 10-minute-distance.

iso10_areas <- iso10_two_modes %>%
  mutate(area = st_area(iso10_two_modes)) %>%
  st_set_geometry(NULL) %>%
  pivot_wider(names_from = mode, values_from = area) 

iso5_areas <- iso_two_modes %>%
  mutate(area = st_area(iso_two_modes)) %>%
  st_set_geometry(NULL) %>%
  pivot_wider(names_from = mode, values_from = area)

Below, I plot the 20 recreation centers by area reachable within 10 minutes. My areas have been calculated in survey feet^2, so I transform my labels into acres by dividing the values by 43,560 (the number of square feet in an acre). By using this number as a factor in my breaks as well, I can plot my gridlines at even intervals.

ggplot(iso10_areas, 
       aes(x = as.numeric(bike), y = as.numeric(walk))) +
  geom_point() +
  geom_smooth(method = "lm",
              linetype = "dashed",
              color = "orange3",
              se = FALSE,
              fullrange = TRUE) +
  scale_x_continuous(name = 
            "\nArea within 10-minute biking distance\nof a recreation center (acres)",
            breaks = breaks <- seq(0, 150000000, by = 6534000),
            labels = prettyNum(breaks / 43560, digits = 4)) + 
  scale_y_continuous(name = 
            "Area within 10-minute walking distance\nof a recreation center (acres)\n",
            breaks = breaks <- seq(0, 10000000, by = (4356000/2)),
            labels = prettyNum(breaks / 43560, digits = 4)) +
  theme_minimal()

And within 5 minutes:

ggplot(iso5_areas, 
       aes(x = as.numeric(bike), y = as.numeric(walk))) +
  geom_point() +
  geom_smooth(method = "lm",
              linetype = "dashed",
              color = "orange2",
              se = FALSE) +
  scale_x_continuous(name = 
            "\nArea within 5-minute biking distance\nof a recreation center (acres)",
            breaks = breaks <- seq(0, 20000000, by = 871200),
            labels = prettyNum(breaks / 43560, digits = 4)) + 
  scale_y_continuous(name = 
            "Area within 5-minute walking distance\nof a recreation center (acres)\n",
            breaks = breaks <- seq(0, 10000000, by = 435600),
            labels = prettyNum(breaks / 43560, digits = 4)) +
  theme_minimal()

We can calculate the factor by which the areas increase when their associated isochrones increase from 5 minutes to ten minutes. To do so, I also filter out the entry for House Park Recreation Center, which was excluded by OTP from the 5-minutes isochrone data.

iso10_sans_house <- iso10_areas[-c(9),]

mean(iso10_sans_house$walk / iso5_areas$walk)
## 5.714484 [1]
mean(iso10_sans_house$bike / iso5_areas$bike)
## 10.1717 [1]

We see that on average, the 10-minute walking isochrones are approximately 5.7x larger than the 5-minute walking isochrones. Interestingly, the 10-minute biking isochrones are about 10.2x larger than the 5-minute biking isochrones. I can visualize these differences in a bar chart. I have to give credit to Julia for the structure of the data frame which goes into the chart below!

iso_areas_2 <- data.frame(rec_center = rep(c("A", "B", "C", "D", "E",
                                             "F", "G", "H", "I", "J",
                                             "K", "L", "M", "N", "O",
                                             "P", "Q", "R", "S")),
                          bikearea5 = as.numeric(iso5_areas$bike),
                          bikearea10 = as.numeric(iso10_sans_house$bike),
                          walkarea5 = as.numeric(iso5_areas$walk),
                          walkarea10 = as.numeric(iso10_sans_house$walk)) %>%
  pivot_longer(cols = c(bikearea5, bikearea10, walkarea5, walkarea10))


ggplot(iso_areas_2, 
       aes(fill = name, x = rec_center, y = value)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_x_discrete(name = "\n Recereation center") +
  scale_fill_viridis_d(name = "\nMode and travel\ntime", 
                       option = "plasma",
                       labels = c("10 min. bike ride", "5 min. bike ride", "10 min. walk", "5 min. walk")) + 
  scale_y_continuous(name = 
            "Area (in acres) accessible by...\n",
            breaks = breaks <- seq(0, 150000000, by = 6534000), 
            labels = prettyNum(breaks / 43560, digits = 4)) +
  theme_minimal() +
  theme()